Finite Temperature Density Matrix Renormalization using an enlarged Hilbert space 
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We apply a generalization of the time-dependent DMRG to study finite temperature properties 
of several quantum spin chains, including the frustrated Ji — J2 model. We discuss several practical 
issues with the method, including use of quantum numbers and finite size effects. We compare with 
transfer-matrix DMRG, finding that both methods produce excellent results. 
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The density matrix renormahzation group (DMRG) 
method provides extremely accurate information 
about the ground state of one dimensional systems. To 
study thermodynamic properties, it was subsequently 
adapted to calculate the transfer-matrix of a ID quan- 
tum system. In the transfer- matrix DMRG (TM-DMRG) 
method 0, the usual DMRG sweeping takes place 
in the imaginary time direction, whereas the thermody- 
namic limit in one spatial direction is automatically ob- 
tained by targeting the maximum transfer matrix eigen- 
value and eigenvector. TM-DMRG gives excellent re- 
sults, but is also technically somewhat more difficult than 
ordinary DMRG, in part because the transfer matrix 
is non-Hermitian. A robust finite-temperature method 
based on the original DMRG method would be very use- 
ful, if only because DMRG and TM-DMRG have slightly 
different strengths and weaknesses. 

In the past year, ideas from the quantum information 
field have greatly extended the capabilities of DMRG 3, 
The first major advance was a new approach to 
real-time evolution within a generalized DMRG frame- 
work (consisting of matrix- product statesQ), and shortly 
thereafter within a standard DMRG frameworkj^ llCll| . 
Subsequently, separately Zwolak and Vidal, and Ver- 
straete, et. al. devised methods allowing finite temper- 
ature DMRG. The approach of Zwolak and Vidal intro- 
duced the idea of a matrix-product description of density 
operators, rather than of wavefunctions. Within the den- 
sity operator formulation, the infinite temperature sys- 
tem is trivial to describe, and imaginary time evolution 
is used to reach finite temperature. Verstraete et al. 
argued that a more efficient procedure is to enlarge the 
Hilbert space with auxiliary sites (called ancillas) , and to 
evolve in imaginary time a pure state within the larger 
space. The auxiliary states act as a perfect heat bath, 
and when traced out give exact thermodynamic averages. 

The ancilla approach is especially convenient from the 
traditional DMRG point of view. A wavefunction de- 
scription is usually more familiar and comfortable than 
a density operator description. The ancillas appear geo- 
metrically as another chain parallel to the first, making 
the system resemble a ladder. More generally, the an- 
cillas form a copy of the original system, doubling the 



size of the lattice. As we discuss below, the ancilla states 
can be given quantum numbers, increasing the efficiency 
of the calculation. In this paper, we apply the ancilla 
approach to study the thermodynamics of several spin 
chains. We consider practical issues such as finite size 
effects, and compare our results with TM-DMRG. 

The use of auxiliary systems to study thermodynamics 
in quantum systems originated as a key idea in thermo 
field dynamics 0, 0, Q| . Let the energy eigenstates of 
the system in question be {n}. Introduce an auxiliary 
set of fictitious states {n} in one-to-one correspondence 
with {n}. Define the unnormalized pure quantum state, 
in an enlarged Hilbert space. 



\m) = E 



(1) 



where fi is the matching state to n, (3 is the inverse tem- 
perature, and |'0(O)) — X^n l'^"') is o^'^ thermal vacuum. 
Then the partition function is 



(2) 



and we can obtain the exact thermodynamic average of 
an operator A (acting only on the real states), as 



{A)^Z{P)-\4,\A\i,). 



(3) 



At (3 — 0, the state -0 is the maximally entangled state 
between the real system and the fictitious system. If 
we change basis from the energy eigenstates n to some 
other arbitrary basis s, ij} is still maximally entangled Ii||, 
= X^s l'^^)- ^ natural basis to use is the site basis, 
where the state of each site i takes on a definite value Si. 
One finds 



1^(0)) =nE I 



(4) 



defining the maximally entangled state of site i with 
its ancilla. 

^From the DMRG point of view, the maximally entan- 
gled state between the left and right blocks would be the 
worst possible state to try to represent: all density ma- 
trix eigenvalues would be equal. If one split the system 
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between the real sites on one side and the anciUa on the 
other, one would have exactly this worst case at /3 = 0. It 
is remarkable that, if one pairs each site with its ancilla, 
and splits the system in two respecting this pairing, the 
infinite temperature state is the best possible state for 
DMRG, with only one density matrix eigenvalue being 
nonzero. This leads to the natural ordering of sites for 
DMRG site-1, ancilla-1, site-2, ancilla-2, etc. Alterna- 
tively, one can group together a site and its ancilla into 
a supersite. Although the dimension of the superblock is 
larger with supersites, next-nearest neighbor interactions 
are not generated by the ancilla, simplifying the time evo- 
lution. In our simulations we have used supersites. 

The essence of the ancilla finite temperature method is 
to start in this local f3 = state, and evolve in imaginary 
time through a succession of temperatures /3. To evolve 
in time, we utilize one of the recently developed time evo- 
lution methods, which perform equally well in imaginary 
time. The most efficient of these utilizes a Suzuki- Trotter 
break-up of the Hamiltonian, and each DMRG step con- 
sists of evolving the state using the link evolution opera- 
tor exp(— ri?i^i_|_i/2) between the two central sites [gllioj. 
This method requires nearest neighbor interactions, at 
least in its simplest form. Alternatively, one can evolve 
in a basis optimized for a single time step by solving ex- 
plicitly the corresponding differential equationjl^], which 
does not require local interactions but is less efficient. 
Note that the Hamiltonian only applies to the real sites; 
the ancillas evolve only through their entanglement. 

In order to utilize conserved quantum numbers, it is 
useful to think of each ancilla as being the anti-site of 
its site. A state of the ancilla is given opposite quantum 
numbers to the corresponding state of the real site. In 
this way, the state of interest has both total charge and 
total z component of spin equal to zero. 

The infinite temperature starting state has a correla- 
tion length of and requires only one state per block. As 
the system evolves in imaginary time, longer range en- 
tanglement is produced and the correlation length grows. 
The number of states needed for a given accuracy grows 
as the temperature decreases. It is most natural to slowly 
increase the size of the basis, in order to keep a roughly 
constant truncation error. One may wish to set a mini- 
mum basis set size to make the early evolution essentially 
exact with little computational cost. In the test calcu- 
lations below we kept the truncation error below 10"^", 
which in the systems considered corresponded typically 
to a maximum of m = 500 DMRG states. It turns out 
that most of these states have = 0, and therefore the 
total size of the basis is of the order of 10^ states for the 
spin-i Heisenberg chain. In more difficult systems one 
would use a less stringent error criterion. 

To illustrate the method, we begin by looking at the 
spin S = 1 Heisenberg chain, using the Trotter time evo- 
lution method. We compare with the TM-DMRG re- 
sults of Xiang li|. The TM-DMRG results have a small. 
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FIG. 1: Specific heat and magnetic susceptibility of the S — 1 
spin chain of length L — 64 obtained with the Suzuki- Trotter 
time-evolution algorithm. We compare with results from TM- 
DMRG in the thermodynamic limit. 



well-controlled Trotter error owing to the formation of 
the transfer matrix; similarly, our time evolution has a 
different small, well-controlled Trotter error. The TM- 
DMRG results are in the thermodynamic limit, whereas 
our results here were on a L = 64 site system with open 
boundary conditions. 

We calculated the specific heat Cv by taking the nu- 
merical derivative of the energy with respect to the tem- 
perature, using energy differences between adjacent time 
steps. In order to avoid edge effects we calculated the lo- 
cal energy in the center of the chain. We also calculated 
the magnetic susceptibility, using the formula 

i 

where the correlations were calculated at equidistant 
points from the center of the open chain, following the 
procedure described in Ref.^^l- We have used half in- 
teger spins at both ends, as in Ref.^^l- Results for 
these thermodynamic quantities are plotted in Fig^ 
The agreement between the our results and Xiang's TM- 
DMRG is very good for both quantities, for temperatures 
down to T w 0.05. At high T we see slight deviations; 
we have checked our results at high T using smaller time 
steps and have found no difference, so we believe the dif- 
ferences are due to Trotter error in the TM-DMRG. 

As a second test example, we choose the spin-i Heisen- 
berg chain with nearest and next nearest neighbor inter- 
actions, with the Hamiltonian 



H — JiSi.Si+i + J2Si.S 



i+2- 



(5) 



Since it is not trivial to use the Suzuki- Trotter break up 
for the frustrated case, we used the time-step targeted 
method[il. 

In Fig[21 we compare our results for Cy and x for 

the unfrustrated chain (J2 = 0) with results from the 
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FIG. 2: Specific heat and magnetic susceptibility of a, S — 1/2 
spin chain of length L — 64, compared to exact L = oo results 
using the Bethe Ansatz. 



Bethe Ansatz calculations of Ref. in the thermody- 
namic limit. The agreement is excellent for the entire 
range of temperatures studied. Finite size effects were 
not apparent down to T « 0.1. 

When frustration is introduced, it is well known that 
this model is gapless for J2 < J2C ~ 0.2411Ji. At this 
value the chain breaks the translational symmetry by 
dimerizing, and an exponentially small gap opens. (20( At 
the point J2 — 0.5 Ji the exact ground states become two 
dimer coverings, and the correlations extend only to one 
lattice spacing. The frustration present in this model 
makes reliable quantum Monte Carlo simulations very 
difficult, due to the appearance of the minus sign prob- 
lem, and the most accurate results for thermodynam- 
ics quantities have been obtained using transfer-matrix 
DMRG. milil 

In Figl^lwe show our results for the specific heat and 
susceptibility for different values of frustration J2, be- 
low and above the critical point J2c- Notice that due 
to dimerization, we have to symmetrize two correlations 
for each distance, {SfS^) and (S'f+iS'|+i). For small val- 
ues of frustration, J2 < J2C, the chain behaves as in the 
unfrustrated case, and the agreement is excellent, even 
for values of frustration up to J2 — 0.5Ji. At this point 
we see more evident finite size effects reflected in the 
curves by a hump that appears in the specific heat at 
small temperatures. A detailed study of this hump shows 
that this is indeed a finite size effect, since it moves to 
lower temperatures as we increase the size of the chain 
as shown in the inset in FigOl At this value of the frus- 
tration, the ground state is degenerate in the thermo- 
dynamic limit, or in chain with periodic boundary con- 
ditions. This degeneracy is lifted in finite chains with 
open boundary conditions. Similar characteristics can 
be observed in Ising chains, where the two Neel configu- 
rations are ground states, when a small off-diagonal cou- 
pling is introduced. 1231 B v imposing a version of smooth 
boundary conditions [24| , where we turn on J2 slowly and 
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FIG. 3: Specific heat and magnetic susceptibility of a frus- 
trated S — 1/2 spin chain of length L — 64 obtained using 
time-step targeting, compared to results from transfer-matrix 
DMRG 113,123. Inset: results for J = 0.5 for different chain 
lengths and smooth boundary conditions. 



smoothly from 0.0 at the edges to 0.5 Ji in the central re- 
gion, we are able to eliminate the hump. 

With our technique detailed spatial correlations func- 
tions are as easily obtained as with ground state DMRG. 
We have calculated the spin-spin correlations and fit 
them to an expression of the form 

C(r) « Aexp(-r/e). (6) 

Figure 0] shows the results for versus T for S* = 1 
and S = 1/2. The agreement with the results from TM- 
DMRG is excellent. For S' = 1 we notice the same min- 
imum observed in the TM-DMRG simulation by Xiang 
16"| . We studied the system using time-targeting and also 
reducing the time step, and we found that the source of 
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FIG. 4: Correlation length as a function of the temperature 
for the 5=1 Heisenberg chain (left panel), and the frus- 
trated Heisenberg chain (right panel), for different values of 
Jil J\- Our calculations are on finite chains of length L — 64. 
We add for comparison results from TM-DMRG in the ther- 
modynamic limit. For 5 = 1 we show results using different 
truncation errors, and at T = from ground state DMRG. 



that minimum can be attributed to the DMRG trunca- 
tion error, as can be seen in the figure. 

To summarize, we have described a DMRG algorithm 
to study strongly correlated quantum systems at finite 
temperatures by using an enlarged Hilbert space with 
ancillary degrees of freedom. We have illustrated its 
application by calculating thermodynamic quantities of 
gapless and gapped systems, including frustration. The 
ideas presented here are simple to implement as an ex- 
tension to standard DMRG codes, and are not restricted 
to nearest-neighbor interactions or to single chains. 
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